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We discuss the response of neutron stars to the tidal interaction in a compact binary system, as 
encoded in the Love number associated with the induced deformation. This problem is of interest 
for gravitational-wave astronomy as there may be a detectable imprint on the signal from the late 
stages of binary coalescence. Previous work has focussed on simple barotropic neutron star models, 
providing an understanding of the role of the stellar compactness and overall density profile. We 
add realism to the discussion by developing the framework required to model stars with varying 
composition and an elastic crust. These effects are not expected to be significant for the next 
generation of detectors but it is nevertheless useful to be able to quantify them. Our results show 
that (perhaps surprisingly) internal stratification has no impact whatsoever on the Love number. 
We also show that crust elasticity provides a (predictably) small correction to existing models. 



I. INTRODUCTION 



When a massive body is exposed to a relatively weak external gravitational field we expect it to respond by changing 

Pm \ shape. This is most easily understood by considering the gravitational effects of the Moon on the Earth. The oceans 

C/^ ■ move to reach equilibrium as the moon orbits, leading to the observed tides. The effect will also deform the Earth's 

r^ I elastic crust, again to reach an equilibrium with the passing body. The deformation of the massive body modifies 

Q_|. the tidal potential, an effect that can be expressed in terms of dimensionless quantities known as the Love numbers. 

1^ ' We will consider this problem for a neutron star in a close binary system, focussing on the tidal Love number; the 

measure of the tidal response of an object to an external gravitational field [l|. 

Specifically, we are interested in ^2, the measure of a tidal deformation by a quadrupole perturbation due to an 
53 \ external gravitational field. This quantity is of direct relevance for future gravitational- wave astronomy, as it provides 
a potentially detectable "finite-size correction" to the late stages of a binary inspiral signal 0]. The effect depends 
on the density distribution of the star, and hence may be used as a diagnostic to use observational data to constrain 
^ , neutron star theory. The implications for observations, and the formalism required to work out the quadrupole 
(^\ ' deformation in relativity, were first discussed in [2|. Since then the formalism has been generalized to determine, 
\^ , in particular, the shape Love number, see for exarnple [3|, \^. Applications have mainly considered barotropic fluid 
^^D ' models, with realistic equations of state studied in [J] and the difference between neutron stars and self-bound quark 
^^ \ stars being quantified in p . 

f — ' The main aim of the present work is to develop the formalism required to study the tidal deformation of real 

^^ [ neutron stars, accounting for key aspects like the elastic crust and variations in the interior composition. These are 

'"^ . important developments, even though we do not expect these features to be easily observed in a binarygravitational- 

wave signal. The Love number leaves an imprint that is borderline detectable by advanced LIGO y. The effect 

will be more important for the third generation Einstein Telescope [6|, but small corrections to it may not seem that 

relevant. Nevertheless, as a matter of principle it is important to move beyond back-of-the-envelope estimates and 

^ , quantify the actual effect. Moreover, developments in this direction are required to address a number of questions of 

j^ ' more immediate astrophysical relevance. By developing the required relativistic perturbation formalism for the static 

deformations induced by the tidal interaction, wc prepare the ground for studies of general crust asymmetries, e.g., 

neutron star "mountains" which lead to gravitational- wave emission at twice the star's spin frequency [7|,|8|. The work 

presented here represents key steps towards modelling such mountains in general relativity, which would allow us to 

consider realistic equations of state for the first time in this context. A closely related problem to that considered here 

concerns the breaking of the crust during binary inspiral. It is easy to estimate (based on the energetics involved [5|) 

that the crust will break before the final merger, but when exactly does this happen and where does the crust yield 

first? Another related problem concerns the quasipcriodic oscillations observed in the tails of magnetar flares [3, [l3l ■ 

Again, a key issue concerns the build-up of stresses in the crust of these strongly magnetized stars, the eventual 

fracture and the induced dynamics. The formalism developed here provides a useful first step for a discussion of this 

problem as well. 



II. NEUTRON STAR RELASTICITY 

The structure of the neutron star crust is relatively well understood [il|, although issues associated with the 
dynamics of the supcrfluid neutron component and the possible pasta phases in the deep crust layers need to be 
explored further. The theory required to model the crust in the framework of general relativity is also well developed. 
Building on pioneering work by Carter [IJ, [ij] and more recent efforts by Karlovini and Samuelsson [l4|, tw o of us 
have recently completed a formalism for Lagrangian perturbations of a relativistic elastic system [1^ (see [13 for the 
corresponding model in Newtonian gravity) . These developments allow us to consider a range of relevant astrophysical 
applications. 

In order to quantify the effect that the neutron star crust has on the tidal deformability of the star, we need to 
formulate and solve the problem for static (quadrupole) deformations of a given stellar model. As usual, this is a 
two-stage process. The first stage is straightforward. If we assume that the unperturbed elastic star is relaxed, i.e., 
that the crust is unstrained, then we simply need to solve the usual Tolman-Oppcnheimcr-Volkoff (TOV) equations 
for a non-rotating relativistic perfect fluid star. The crust manifests itself only at the linear perturbation level. At 
this second stage, we need to consider the static perturbations of a model that accounts for the associated elastic 
stresses. 

In this section, we develop the formalism required to solve the tidal deformation problem for elastic relativistic stars. 
In particular, we develop the relevant perturbation equations and discuss the key differences from the perfect fluid 
case. Throughout most of the analysis, we use geometric units where c = G = 1, and all primes denote differentiation 
with respect to the radial coordinate. We restore the physical units only when they are useful for analysis purposes. 

A. The background problem 

The equilibrium of a static, spherically symmetric, relativistic star is described by a spacetime metric gab given by 

The fluid four- velocity is simply given by 

u° = e-^/'^t" , (2) 

where t° = [dtY is the timelike Killing vector of the spacetime. 

We want to model a neutron star with a fluid core, an elastic crust and possibly a fluid ocean. The corresponding 
equilibrium problem simplifies significantly if we assume that the crust is relaxed in the unperturbed configuration. If 
this is the case, then the problem reduces to that of a fluid system. Physically, this assumption makes sense provided 
the crust has had time to release any built up strain, e.g. via plastic flow, before the tidal interaction with the 
star's companion becomes sizeable. The upshot is that we can use the perfect fluid stress-energy tensor to model the 
background star. In other words, we have 

Tab = (P + P)UaUb + Pgab = PUaUb + P -Lab , (3) 

where _Lq6= gab + UaUb is the usual projection orthogonal to the fluid frame. In the particular case of a barotropic 
star, the pressure and energy density are related by an equation of state of form P ~ Pip)- We will discuss the more 
general case, where the (cold) equation of state is such that pressure also depends on the composition of matter at 
supranuclear densities, later. 

Once the equation of state is provided we need to solve the standard TOV equations starting with 

X' = -^{M~47Tr'p) , (4) 

where the gravitational mass Af (r), inside radius r, is obtained from 

M' = 47rrV . (5) 

Then using 

A o 



we arrive at 



:^(Af + 47rr3p) = -^-^P', (6) 

^ (p + P)(M + 47rr3p) 

r{r - 2AI) ' ^ ' 



These relations provide enough information to determine the background stellar model. 



B. Static perturbations of a fluid star 



In order to investigate the tidal deformations, we need to consider static perturbations of the system. Since we plan 
to account for elastic contributions, it is natural to approach the problem within Lagrangian perturbation theory. The 
relevant formalism has recently been developed J15l | . For the present purposes it is worth noting that the problem we 
consider is very similar to that of the zero-frequency subspace of perturbations for slowly rotating stars. This means 
that we can draw on work relating to the gravitational-wave driven r-mode instability. In particular, the spacetime 
part of the perturbation equations is identical to that considered by [17| . 

Due to the nature of the tidal interaction, we focus on the polar perturbations. These lead to the electric-type Love 
numbers, as discussed in |3,I3- We use the standard Regge- Wheeler gauge, representing the perturbed metric hat by 
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(noting a typo in the off-diagonal element in [i7[)- Due to the spherical symmetry of the problem, the different 
multipoles do not couple. In fact, we can, without loss of generality, assume that the perturbations arc axisymmctric 
(set TO = 0). 
Wc introduce a static displacement vector given by |3l| . 



^« = -W{r)Yi''\'' + V{r)V''Y{' 
Within rclativistic Lagrangian perturbation theory [IJ, [l^, |2C)| - (23 | we then have the perturbed four- velocity 

where 

Apab = hab + 2V(aCb) ■ 

Moreover, the Lagrangian perturbations arc related to the Eulerian perturbations by 

A = (5 + £5 , 
where C^ is the Lie derivative along ^''. It follows that the perturbed 4- velocity, (5u", is given by 



(9) 

(10) 

(11) 
(12) 

(13) 



Since the displacement vector is taken to be static, we have 



5u* 



-e-^'/'^HoY^"' , and 5u^ = 



(14) 



where j = 1 — 3 represents a spatial index. 

Wc also know that the perturbed number density can be obtained from the displacement vector, since [22 



An = 5n + ^V aU 



iyi ^'^^ A,9,b . 



Explicitly, we have 



n 



K + l-H-A - l{l + l)V + rW' +il + ]-r\' ] W 



Yr 



(15) 



(16) 



In the case of barotropic matter the energy density will be a function only of the number density of the constituent 
particles, p — p{n) which means that. 



Ap ~ ——An — fj,An , 
an 



(17) 



where we have introduced the chemical potential fj, [22|. For algebraic simplicity later we define 

2 ' 



-^9- 9 



r^ ( K + I-H2] - l{l + l)V + rW + ( 1 + \r\' ] W 



(18) 



Combining (|16p with the Gibbs relation P + p = nfi, we sec that 



Ap^--iP + p)±,Yr. (19) 



Finally, in the case of a barotrope, the perturbed pressure follows from 

rip fip 

AP = -—Ap = clAp , and 6P ^ -—5p , (20) 

dp dp 

where we have defined the speed of sound, c^. These relations show that we can choose to work either with the 
components of the displacement vector or the perturbed density/pressure. One of the variables, (5p, W and V is 
redundant. The situation is, of course, different for a non-barotropic model. We will discuss this case later. 

To complete the specification of the fluid problem we need the perturbed Einstein equations (for I > 2). The 
right-hand side of the equations follows from ([3|) . In the case of a perfect fluid we have 

STJ' = {5p + 5P)UaU^ + 5P6^ + {P + p) {u''6ua + UaSu'') , (21) 

where the scalar perturbations are to be expanded in spherical harmonics. It is worth noting that all off-diagonal 
components vanish identically since Su^ = 0. 

The corresponding left-hand side of the Einstein equations, expressed in terms of the perturbed metric, obviously 
does not depend on the matter content of the spacetime at all. Basically, the required expressions can be found in a 
number of other studies, including |17l. Il8| and |23j . In the case of a spherically symmetric system, we end up with 
six coupled ordinary differential equations to solve. The final fluid equations are easily obtained from the full elastic 
equations discussed below, so we will not give them here. 

C. Elasticity 

Having outlined the fluid perturbation problem, and the foundations of Lagrangian perturbation theory in the 
relativistic setting, we can move on to consider the role of the neutron star's elastic crust. As we are working with 
the Einstein equations, we only need to consider the relevant alterations to the stress-energy tensor. A perfect fluid 
cannot support shear stresses and hence does not introduce off-diagonal stress-energy terms. When we account for 
elasticity we need to include such stresses [15|, |la| . 

As discussed above, we assume the background star is in a relaxed unstrained state. Given this assumption, the 
elastic crust docs not contribute to the stress-energy tensor of the background star, which leaves our equilibrium 
equations unaltered. The elastic contributions enter only through the perturbed stress-energy tensor. Following [l5|, 
we have the Lagrangian perturbation of the anistropic stress tensor; 

ATTab -= -2pASab , (22) 

where p is the shear modulus (not to be confused with the chemical potential), and 

2Asab = U^l-^b ~l ^abl-'A Ag,d ■ (23) 

In the case of an unstrained background, these expressions lead to 

S^a' = -2A (l-^l-'' -I K'l-'A Affcrf , (24) 



which should be added to the fluid result (|21l) . For future reference, it is useful to note that 

dTTt'' = . (25) 



In order to complete the perturbed Einstein equations in the elastic case, we need to express the perturbed stress 
tensor in terms of our chosen variables. Thus, we find that 



S-K^ 



-r^jl [r^{K - H2) - 1(1 + 1)V - 2rW' + (4 - r\') W] 



(26) 



(here, and in the following, we suppress the angular dependence for clarity). Moreover, since the anisotropic stress 
must be trace- free, c.f., p4| . it follows immediately that 



S-Ka 



Stt^^ 



-5nJ . 



Note that this result requires that the background is relaxed. That is, we have 



Stto 



-^< 



-fi [r^(H2 - K) + l{l + l)V + 2rW' - (4 - rA') W] 
o 



We also find that 



and, finally. 



Stto" - Sn^^ 



SfiV , 



Stt^ 



2/i 



{rV' -2V + e^W) 



(27) 
(28) 
(29) 
(30) 



Given the above results, the perturbation equations for the stellar interior become more complicated. However, the 
[t <]-component of the perturbed Einstein equations is unaffected by the presence of the crust, so we have from the 
fluid case, 



e-\''K" 



-M 3 _ ^:^] ,K' 



l{l + l)-l 



e-^rH'^ 



l{l + l)+e-^{l-r}^) 



H2 = -Sirr^Sp . (31) 



Moreover, because our star is spherically symmetric and (57r/' = 0, we find that Hi — (c.f., (23). Finally, because 
Sit J' is traceless, we can take the trace of the perturbation equations to arrive at a second equation without explicit 
dependence on the elasticity. This equation takes the form 

- r'^Hi; + (-r\' - rv - 2 j rH'^ + r^v'K' - -r'^v'H'^ + l{l + l)e^Ho 

+ [2 (e^ - 1) - r (A' + 3i^')] ^2 = Snr'^e^iSSP + 5p) . (32) 

In practice, it is useful to introduce two additional variables closely related to the traction. These variables, which 
vanish in the fluid case, are key to the implementation of the relevant junction conditions at the crust-core interface. 
Hence, we define a variable linked to the radial component of the traction; 



and its 9 component; 



Ti = -h [r''{K - H2) - l{l + l)V - 2rW' + (4 - rA') W] , 



T, 



(33) 



(34) 



^4/1 {rV - 2V + e^W) . 

Working out the perturbed Einstein equations, following the same strategy as in the fluid problem (see [l7[ and 
|23|). we find that the difference between the [6 9] and the [(p (p] components leads to 



H2- Ha = M-KJiV 



Meanwhile, the [r 9] equation leads to 



K' 



-. — I' r^f 



e-iJo]' + 5^(r^' + 2)F-8^^ 



r r 

In terms of the traction variables we can write the sum of the [9 9] and the [ip ip] component as 



SP^ 



1 



-,-A /,,/ 



,-A 



16nr 



{v' + A')-ffo + ^r {2^6^ [2 - l{l + l)]V- - W - A') + 1 



e^ % T2 

T2 + —Ti - + — 

2 2 2 



(35) 



(36) 



(37) 



where we have also used (|35)) and (|36)) . 

Finally, the [r r] component of the Einstein equations leads to 

[l{l + 1) - 2] e^K = r^v'H'^ + [l{l + l)e^ - 2 + r^v'f - r {v' + A')] H^ 

+ 327r/i p (i^')^ + /(/ + 1) - 2J y - 47rr (;/' + A') T2 + 87r(rT^ - T2) - 247rTi . (38) 



It is worth noting that, by adding (|37p and (p8| then subtracting (|3ip we have the trace of the Einstein equations, 
i.e., we recover (l32l). 



D. The final perturbation equations 

The set of equations that we need to solve in the crust region combines the conservation law (|16p with the perturbed 
Einstein equations (|35p ~ ((55)) . including the definitions ([55]) and ([M]) . We want to formulate the problem in such a 
way that the required integration becomes as straightforward as possible. To do this, we note that we can reduce the 
order of the system. However, in the elastic case this reduction does not lead to the same level of simplification as in 
the fluid problem, where we only need to solve a single second order equation for iJo, C-f-, dM]) below. This is of course 
as expected, as the problem now has additional degrees of freedom due to the elasticity. For example, it is natural 
that we will have to solve for the components of the displacement vector. We take the view that these components 
are obtained from the definitions (|33|) and (p4| . That is, we integrate 



W'-[\-^^W=\{K-iI,)- 



32TTJlr 



Ijl + l) 
2r 



^-8^^- 



and 



r r 4:fl r 



(39) 



(40) 



Next we note that we have several relations involving the perturbed pressure. If we focus on the barotropic problem 
it follows from (TTSl) that we should have 



where 



AP = c',Ap = -^iP + p)cl±g, 



W 1 

AP = SP + —P' = SP (P + pWW . 

r 2r 



(41) 



(42) 



Using ((39)) to remove W' from _Lg we arrive at the algebraic relation 



SP = -^^c^\lr-K-hl + l)V+U-'^]w^lT, 



2c 



8/i 



We can rearrange this as an equation that determines Ti , giving 



^ 8/i 



{P + p)c 



^SP+-r^K--l{l + l)V+[i-—^]W 



(43) 



(44) 



By combining (J37]) and psp . in such a way that Tj is removed, we arrive at a second algebraic relation. This 
provides another expression for (5P, giving 

IQirr'^e^SP = r'^v'H'^ + [2 - l{l + l)]e^K + \l{l + l)e^ - 2 + r^ [v'f] Hq 

+ 32TTJ1 1?'^ {v'f + [l{l + l)-2]{l-e^)\v - Sn [rv' + 2) T2 + 87r (e/ - 3) Ti . (45) 

We can use this for the left-hand side of ([37)) . which then becomes an equation that we can integrate to get T2; 



T' 



5<-'-^') + ; 



1 Au c^ 

Ta = -2eV(5P + tt i^' + A') Hq + -^e^ [2 - l{l + l)]V+ —Ti . 
STT r r 



(46) 



Finally, we have (|36p for K' ^ which we can use in (|32)) (together with ((35|) for H2) to get 
1 



r'H'^ 



^r{\'-v')-2 



tH'q + {/(? + l)e^ + 2(e^ - 1) - r(A' + 3;^') + r^ (z^')'} ^0 



87r<^r^e-^(35P + (5p) + 16A 



l-e^ + rf^' + ivj-iKf 



y + 4r2i/' (/iy)' + ri^'T2 I . (47) 



These equations completely specify the perturbation problem in the elastic region. It is easy to verify that we have 
the same number of equations as we have unknowns. 

E. Interface conditions 

Having determined the perturbation equations for the elastic crust region, we need to connect them to the fluid 
perturbations for both the core and the fluid ocean near the star's surface. This requires a set of interface conditions 
at a radius r^, which could represent either the crust-core or the crust-ocean interface. The interface problem has 
already been analyzed in detail, the most relevant work being that of Finn [l9[. As our analysis is identical, we simply 
summarize the results here. 

From our assumption that the unperturbed crust is in a relaxed state we know that all background quantities are 
continuous across a fluid-crust interface. This includes the density. It may, of course, be that the true equation of 
state is such that the interface is associated with a small density discontinuity (c.f., the discussion in Q). We do not 
account for this possibility here, but it would be very easy to do so, should it be required. 

To determine the behaviour of the perturbed quantities across the fluid-elastic interface we calculate the intrinsic 
curvature at the interface. From the Israel junction conditions we know that the intrinsic curvature must be continuous. 
In the fluid problem (in fact, even for multi-fluids j24l|). this results in the continuity of all perturbed metric variables, 
as well as their first derivatives. The proof of this relies on the fact that H2 = Hq in the fluid case. This is, however, 
no longer true when we consider the elastic problem. Consequently, the behaviour at a fluid-elastic interface is a little 
bit more complicated. 

Nevertheless, we know from the analysis of Finn [l9[ that the continuity of the first fundamental form demands the 
continuity of ly, Hq, Hi and K at the interface. The first of these conditions follows from the continuity of ^'' = W/r, 
which tells us that in order to "avoid a void" , the radial displacement component W must be continuous across the 
interface. Note also that the condition on Hi is irrelevant to us as this variable vanishes identically in the problem 
under consideration. 

Combining (j35p with the stated continuity conditions we see that we should have 

[H2l^^?,2i:[jlVl^ , (48) 

where [A\r^ is shorthand for lim£^o^(^c + e) — liin£^o^(^c — e)- This result shows that we should expect a jump 
in the H2 perturbation; first of all the shear modulus will vanish sharply as the crust gives way to the fiuid core, 
and secondly there is no reason why the tangential displacement V should be continuous (at least not as long as we 
are ignoring viscosity and magnetic field stresses). In the limit where /i — > we obviously recover the standard fluid 
interface condition. 

The remaining interface conditions can be obtained either by imposing the continuity of the extrinsic curvature 
Kab^ or the surface stresses Sab on the perturbed hypersurface. The two approaches are related since |25| . 

Sab = ^{[Kab]-[K]^''^gab) , (49) 

where '^'gafc — dab ~ NaNb is the induced three-metric on the surface and K — Kab Vab, and Na is the normal to 
the surface. Since the right-hand side must be continuous across the surface we demand that [Sab] ~ 0. This provides 
two additional conditions [19|; 

[Ti+r2AP]^^=0, (50) 

for the (perturbed) radial traction, and 

[^2],^ = , (51) 

for the horizontal part. Since the background variables are taken to be continuous across the interface, and since the 
intrinsic curvature conditions lead to [W]r^ = 0, condition ([50]) yields 

[5P + Til =0. (52) 



We now have six conditions at each fluid-crust interface, to combine with the six ODEs in the crust. The problem 
is therefore well posed. 

The implementation of the interface condition depends, to some extent, on the physical model considered. We are 
interested in a stellar model with a fluid core and an elastic crust that transitions back to a fluid ocean near the 
surface. Consequently, we need to impose the junction conditions at two fluid-elastic interfaces. Integrating from the 
centre of the star, we solve the fluid problem up to the crust-core interface. At this point the conditions provide us 
with all information needed to continue the integration using the elastic equations. Finally, the same conditions are 
imposed at the crust melting point. Here, the continuity of Hq and K implies that we have the information required 
to integrate the fluid equations to the actual surface of the star. 

Specifically, we first of all have the continuity of Hq and K at each interface. Moreover, from the condition [T2] = 
we realize that since the shear modulus p, has a finite value in the crust, but vanishes in the fluid, we must have 
T2 = at the interface. The condition [W] = is a little more complicated. Since we do not calculate the radial 
perturbation, W, in the fluid we need a means to initiate the integration in the crust. To do this, we make use of the 
condition ([50]) . Considering the algebraic relations between SP and Ti, Eqns. (|43| and (|44|). and using continuity we 
have 



le-rrr'^e^SPE = r^v'H'j^ + [2 - l{l + 1)] e^Ke + [l{l + l)e^ - 2 + r^{v'f] He 

+ 327rA {r^{v'f + W + 1) - 2](1 - e^)} Ve + 87r(e^ - Z)Tie , 



(53) 



where the subscript E denotes a quantity calculated in the crust region. Using the jump conditions we can express 
He and Ke in terms of the fluid quantities. We also use 
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-r^K 



e--1{1 + \)Ve 



TV , 

^-^)We 



(54) 



The continuity of the radial traction closes this set of equations: since Tip = (F indicating a variable in the fluid 
region) identically we must have 



Tu 



{p + P)Hf - SPe 



(55) 



Combining equations 
base of the crust. 



53 |) - (|55p we have the information required to determine We, the radial perturbation at the 



III. THE LOVE NUMBER 

In order to implement the formalism in an astrophysically meaningful context, we will quantify how the crust 
elasticity affects the tidal deformation of a neutron star. This effect can be expressed in terms of the tidal Love 
number. A static spherically symmetric star of mass M and radius R exposed to a time- independent external tidal 
field £ij will develop a quadrupole moment Qij. To linear order we relate this quadrupole moment Qij to the tidal 
moment £ij thus defining the Love number, ^2, [23|, 



Q^ 



:k2R''£^ 



(56) 



We briefly review the procedure used to calculate the tidal Love number below. For a detailed description we refer 
the reader to any of 0-11] or [H]. 

Following [2^1 , the Love number is extracted from the asymptotic behaviour of the gravitational field of a tidally 
deformed body. To do this we note that the vacuum perturbation problem reduces to a single ODE for Hq 



H'^ 



=^,-X']H'q 



l{l + l)e^ 



(A' 



/\2 



Ho = 



(57) 



where we have used v = —A and M{r) = M = const, exterior to the star. This equation may be solved in terms of 
the associated Legendre polynomials [2, 0, [i3j leading to; 



H{r) = apPa (]^ - l) + ^qQn (^ - 1 
Wc have recorded both the decreasing solution Pim{x) and growing solution Qim{x) for to = 2. 



(58) 



Since the problem is studied within linear perturbation theory the amplitude of the solution is arbitrary, so we 
demand that the function 



y{r) ^ r 



H'jr) 
H{r)- 



matches across the surface at the star. Substituting the general solution Eqn. ([55]) into Eqn. ([CT)) . we get 



y{x) = (1 + x) 



Pi2{x) +aiQi2{x) 



(59) 



(60) 



where x = R/M — 1. Following, [j| we have defined a; = aq/ap which is determined by matching Eqn. (j60p across 
the surface of the star. This leads to 



ai 



Pl^{x) - CyiPnix) 



(61) 



R 



Q'i^{x) - CyiQnix) 
where C = M / R is the compactness of the star. Substituting I — 2, and taking the asymptotic expansion we get 

^o«3(-) aP + 3(-) «,, (62) 

which is the result used in j23j . 

Finally, we relate the coefficient a; to the Love number by comparing to the response of a spherically symmetric 
star to an external quadrupolar field, £ij [23l.l26lj. 



1 — gtt M 3Qij fx^x^ 1 



2 r 2r3 V ^ 3 

where we have dropped terms of order 0{l/i'^) and 0{i'^). Using 



1 



-S'^ ] + -e^.x'x^ , 



9tt 



1 - ^"j (1 - H^Yi^) 



and Eq. (|56p . we determine the Love number, ^2, to be 



46 /A/\'' 



(63) 



(64) 



(65) 



where ai is determined by Eqn. (j6ip from the solution to the interior problem. 

In the perfect fluid problem, i.e., ignoring the elastic contributions, the interior equations can be reduced to a single 
second order equation for Rq [2j|, 



r^W' 



h>'-''^ 



r'H' 



2 (1 - e^) - l{l + l)e^ + 2r{2iy' + A') - r^ [v'Y Ho = -8Trr^e^{dP + Sp) . (66) 



Combining this equation with 6P = c^Sp (for a barotropic model) and (|37p . obviously still ignoring the elastic 
contributions, we have a complete formulation of the problem, and a route to determining the coefficients ap and ag 
required in the exterior. When considering the elastic system, we have to solve the set of equations discussed in the 
previous section in order to determine ap and ag. 

For future reference, it is worth noting that (|66p contains the same information as p2p. as is necessary for the 
problem not to be overdetermined. However, one can show that this implies that 



Snr-^e^i^'Sp 



r^X" - r^ (A') + 2 (1 - e^) 



i^o- 



(67) 



Making use of the background equations and the sum of the [6 6] and [ip ip] components of the Einstein equations we 
find that we must have 



p'SP = P'5p 



(68) 
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This is (obviously) true for barotropcs, where P = P{p), so the problem is well posed. 

Finally, using Eqn. (|59p and the compactness parameter C, we determine the coefficient 02 and obtain the Love 
number from 

2^{l-2Cf[2 + 2C{y-l)-y]jx 

{2C(6 -3y + 3C{5y - 8)) + 4C^[U - lly + C{3y - 2) + 2(7^(1 + y)] 

+3(l-2C)2[2-y + 2C(y-l)]log(l-2C)}"' . (69) 

IV. RESULTS 

The previous sections specify the calculations needed to account for the crust elasticity, determine tidal response 
and quantify the Love number fc2 for realistic neutron star models. We will now show how the detailed interior physics 
affects the result. We consider both composition variations and the role of the elastic crust. 

A. Stratification 

It is well-known that variations in the composition of the neutron star core, e.g. represented by a varying proton 
fraction, may have repercussions for the global dynamics. In particular, the associated stratification may lead to the 



presence of g-modes in the star's oscillation spectrum j27| . From a mathematical point of view, these effects arise from 
the equation of state depending on two (or more) parameters. The impact on the dynamics depends on the detailed 
timescales of the problem. We will consider two extreme limits. In the first case, when reactions are very fast, the 
perturbed fluid elements will adjust to their surroundings very quickly. In other words, they lose their original identity 
as the system evolves. In this limit one would expect the problem to remain effectively barotropic. In the opposite 
limit, reactions act slowly compared to the dynamics. In the present context, this should be taken as meaning that 
the reactions are slow compared to the binary inspiral timescale. If we focus on the late stages of evolution, which are 
key from the gravitational-wave observation perspective, then this timescale would be of the order of a few minutes. 
If the involved reactions are slower than this, then a perturbed fluid element must retain its identity (composition) 
through the evolution. We expect this limit to be relevant for binary neutron stars, c.f. j27| . 

Let us first analyze the case where the relevant reactions are much faster than the tidal dynamics. In this case it 
is natural to consider an equation of state such that the pressure depends on two parameters; the energy density p 
and a parameter 13 that represents the deviation from chemical equilibrium. Since the fluid elements will have time 
to cquilibriate with the surrounding fluid, they will remain at (local) chemical equilibrium and we will have A/3 — 0. 
However, /3 = also in the background, which means that 

<5/3 = A/3 - rVa/5 = , (70) 

as well. The upshot of this is that the system behaves as a barotrope. The fast interactions do not affect the tidal 
perturbations, and cannot have an impact on the Love number. 

In the opposite case, when reactions are very slow, it is natural to consider an equation of state with two charge 
neutral components, neutrons (with index n) and a conglomerate of protons and electrons (with index p), such that 
P = P{n, Xp) where n = rin -I- rip is the total baryon number density and Xp = Up/n is the proton fraction. Assuming 
that each species is conserved, i.e., that the relevant reactions are too slow to equilibrate the matter, but the two 
components are locked together (i.e. there is only one displacement vector [la]), we replace ([TS]) by 

An^ = --72^ l"'' Agab , x = p, n . (71) 

Using these relations, we find that An is still given by ([T5|) and we also see that Axp = 0. These results simply 
represent the assumption of frozen composition. It also follows that 

AP = ( ^ ) An , and Ap = f |^ ) An . (72) 



dn J ^ \dn 

That is, we have 



^P-[%) Ap. (73) 
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In other words, the Eulerian variations are related by 



SP^ 



dP 
dp 



Sp + AaC 



(74) 



where 



Aa = S/aP = 



dp 
dp 



"^aP 



(75) 



where /3 has the same meaning as in the case of fast reactions. The magnitude of A" defines the Schwarzschild 
discriminant j27| . 

Let us now consider the perturbed Einstein equations in the case of frozen composition. Since we made no as- 
sumptions about 6P and Sp when we discussed the perturbation equations, the analysis in Section [II Bl must remain 
unchanged. In particular, the final differential equation for Hq remains unaltered. However, the fact that this equation 
can be derived in two different ways is now important. As we have already mentioned, the problem is overdetermined 
unless (|68p holds. Combining this condition with ([71|) we see that we must have 



dP 



dP 

dp/ p 



Ap = 0, 



Ap = 



(76) 



That is, the Lagrangian variation of the density must vanish identically. Moreover, we are led to the (possibly 
surprising) conclusion that internal stratification has no effect on the Love number. 

These arguments suggest that Love number measurements will only reveal limited information about the equation 
of state. We may be able to constrain the supranuclear equation of state, e.g. in terms of the inferred compactness, 
but not probe the detailed composition. 



Numerical results: The crust 



In order to quantify the role of the crust elasticity, we need to study the problem numerically for specific neutron 
stars models. In this first proof-of-principle study we will consider a simple model, based on a polytropic equation 
of state and an analytic model for the crust's shear modulus. In principle, we are ready to implement more realistic 
descriptions of both the bulk behaviour and the elasticity, but before doing this it makes sense to develop the 
computational strategy for this simpler model problem. An obvious advantage of the simplified description is that it 
is relatively straightforward to write down a set of non-dimensionalised perturbation equations, see Appendices lA II 
and EH 

In order to facilitate a direct comparison to previous work 0, |j, l2j] we will assume that matter is described by the 
polytropic equation of state. 



P^Kp 



l + l/JV 



(77) 



where K and N are constant. In addition, we need to specify the shear modulus fl. The outer regions of an 
astrophysical neutron star are composed of i) a thin fluid ocean, below the density at which the crust melts for the 



given temperature (typically, p < 10^ g cm 



^), ii) the outer crust, reaching up to neutron drip (10^ g cm 



10^^ g cm ^), and iii) the inner crust, where superfluid neutrons permeate the nuclear lattice (10 



11 — s 

g cm ■^ 



< P< 

< p < 

10^^ g cm '^) Moreover, the bottom layers of the crust may be in the so-called pasta phase with rather different elastic 
properties [Il|, ^^. We will not consider this possibility here. Our simplified model has a single elastic crust of 
thickness ARc. We will present results for a set of (moderately realistic) polytropic models for which the crust region 
starts at 2 X 10^"' g/cm , and stops at 10^ g/cm . The non-crust regions are treated as perfect fluids. In the crust we 
implement a simple linear shear modulus |ll| . 



p = kP + p.a, 



(78) 



where k is a scaling constant and /ip is a constant allowing us to consider an adjustable shear modulus at the top of the 
crust. Both K and ^o are tuneable parameters. The equations presented in Sec. Ill Dl are, of course, independent of the 
specific form of the shear modulus (|78p . Hence, the implementation of more realistic models would be straightforward. 
We have, first of all, tested our numerical implementation by comparing the results for a fluid model to those of 
[3| and [23| , achieving good agreement and hence confirming that the Love number is smaller for larger values of the 
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1,4 




Stellar Radius (km) 

FIG. 1: The mass-radius curve for the stellar models considered in Table[l] (represented by triangles), showing that the considered 
models are stable to radial perturbations. In comparing to the data in the Table it is useful to recall that IMq = 1.4773 km 
in geometric units. 



Po (km-2) 


C 


^2 crust 


^9. fliiiH 


^k2/k^^ fl.iiH 


Ai?c (km) 


R (km) 


M (km) 


3.378893O-03 


0.24813001129(7) 


0.02413726613(6) 


0.024138574(6) 


-5.420(7)e-05 


0.4412458028(3) 


7.9838719798(7) 


1.9810382445(6) 


3.044717C-03 


0.23938571005(1) 


0.02748474657(2) 


0.027486580(6) 


-6.672(4)e-05 


0.4847246944(3) 


8.2035091998(7) 


1.9638028747(2) 


2.710540C-03 


0.22915622628(4) 


0.03172353529(9) 


0.031726271(2) 


-8.623(6)g-05 


0.5398838957(1) 


8.4494607108(1) 


1.9362465306(3) 


2.376364C-03 


0.21709115810(1) 


0.03719492461(6) 


0.037199305(6) 


^1.177(7)e-04 


0.6117312291(5) 


8.7266756593(2) 


1.8944841252(5) 


2.042188C-03 


0.20272655170(3) 


0.04441587924(3) 


0.044423470(8) 


-1.70S(9)e-04 


0.7085096707(6) 


9.0414350438(4) 


1.8329389488(8) 


1.708012O-03 


0.18543653698(4) 


0.05419576631(7) 


0.054210160(9) 


-2.655(3)e-04 


0.8447708420(0) 


9.4018321596(9) 


1.7434431969(9) 


1.373836C-03 


0.16435873698(6) 


0.06785042152(8) 


0.067880961(4) 


^4.499(0)c-04 


1.0487697824(5) 


9.8184762317(9) 


1.6137523525(9) 


1.039659C-03 


0.13827570668(3) 


0.08761155990(9) 


0.087687716(3) 


-8.684(9)e-04 


1.3835220344(5) 


10.305549472(1) 


1.4250071360(2) 


7.054831C-04 


0.10541985896(6) 


0.11742663334(5) 


0.117676862(9) 


-2.126(4)e-03 


2.0250735506(0) 


10.882441793(5) 


1.1472254790(8) 


3.713069C-04 


0.06313942900(5) 


0.16404412956(2) 


0.165596266(2) 


-9.373(0)0-03 


3.7551538607(9) 


11.576364548(1) 


0.7309250475(2) 



TABLE I: A sample of numerical results comparing the Love number for fluid models and the elastic crust models developed 
in this work. The results are obtained using an A'' = 1 polytrope with K = 100 km^. The stellar models are determined by 
the central density pc- We provide the resulting compactness, C — M/ R, mass, M, radius, R, and crustal thickness, A_Rc, for 
the background star. The tidal Love number for both the crust, fc2 crust, and purely fluid models, fc2 fluid, are shown. From the 
differences between the final Love numbers, Afc2 = fc2crust — fc2fluid, we see that the crust produces a very small correction to 
the tidal Love number. For this table we used the elastic parameters k — 0.015 and /xq = 10~^* km~^. 



polytropic index. The tidal response is weaker for stars that are more centrally condensed, which is natural. Adding 
the crust to these models, we obtain the results given in Table HI The results were obtained for an TV = 1 polytrope 
with K = 100 km^ and the crust model ([75)1 . The numerical results given in TableUwere obtained for a shear modulus 
K = 0.015 and a shear constant /xq = 3 x 10~^* km^ , but we have confirmed the general behaviour by varying these 
parameters. The chosen models are all stable to radial perturbations, c.f., the mass-radius curve in Figure [TJ As 
expected, the presence of the crust has only a small impact on the tidal Love number. By comparing the relative 
change in the Love number with increasing central density to the corresponding change in crust thickness (the left 
and right panels of Figure [21 respectively), we see that the effect increases as more of the star becomes elastic. This 
is, of course, as expected. 

The main message is that we have successfully implemented the crust, and the numerical solution is accurate enough 
to distinguish the small effect that the elasticity has in this problem. The qualitative behaviour of the results is easily 
explained. The results in Table |I] show that presence of the crust decreases the tidal Love number. This is expected, 
since the Love number is a measure of the response to the presence of an external gravitational field. In a fluid model 
we expect a larger distortion since fluids do not sustain shear stresses. Crustal models are able to resist deformations, 
and arc thus expected to have smaller Love numbers. 
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pJ{W g/crri») 




0,0001 =- 



pJ{W g/cra') 



FIG. 2: Left panel: The ratio of crust thickness to stellar radius as function of the central density. The results show that 
the crust occupies a much larger fraction of lower density stars. Right panel: The relative change in Love number due to the 
presence of the crust as function of the central density. As the central density approaches the crust-core transition density, 
2 X 10^* g/cm , the contribution to the Love number becomes larger due to the fact that the crust occupies a much larger 
fraction of the star. Triangles correspond to the models listed in Table [H 



CONCLUSIONS 



We have extended the discussion of tidally deformed relativistic stars by including the effects of internal composition 
stratification and the presence of an elastic crust, relevant for mature neutron stars. The most important development 
concerns the formalism required to account for more realistic neutron star physics. Building on a recent extension 
[15J of relativistic Lagrangian perturbation theory to the multi-fluid setting (allowing for one of the components to 
be elastic), we have formulated the tidal deformation problem for realistic neutron star models. The final model is 
obviously still "incomplete" , as we are ignoring the magnetic field and we are not accounting for superfluidity (which 
should be relevant both in the crust and the fluid core), but it nevertheless allows us to consider astrophysically 
relevant questions. 

This paper should be seen as a first, proof-of-principle, study of the relevant issues. We have shown that (perhaps 
unexpectedly) the tidal deformations arc not affected at all by composition variations in the star's core. Having 
implemented the formalism required to account for the crust elasticity, and solved the problem numerically for a 
simple model problem (based on a polytropic equation of state and a simplistic model for the crust shear modulus), 
we have also demonstrated that the presence of the crust has a (predictably) small effect on the tidal Love number. We 
have considered how this effect varies with the crust parameters, and confirmed that the results agree with intuition. 

Prom a formal point-of-view, these are important developments, but it should be stressed that we do not expect the 
new features to leave an observable imprint on a binary gravitational- wave signal. This is fairly obvious since the Love 
number leaves an imprint that is barely detectable by future generations of gravitational-wave detectors in the first 
place y and [6|, and we are quantifying "small" corrections to it. Having said that, it is clear that we need to move 
beyond back-of-the-envelope estimates in this problem area and develop the computational technology required to 
model real neutron stars and actual astrophysical scenarios. The present work represents an important step towards 
this goal. By developing the required relativistic perturbation formalism for the tidal interaction problem, we lay 
the foundation for work on general crust asymmetries, e.g., neutron star "mountains" relevant for gravitational- wave 
astronomy. This problem has not yet been considered in general relativity, which is required if we want to use realistic 
equations of state. A closely related problem concerns the fracture of the crust during binary inspiral, the effect that 
this may have on a detected gravitational-wave signal and possible associated electromagnetic simatures. Another 
important problem concerns the quasiperiodic oscillations observed in the tails of magnetar flares [3, llfl l2Ml • I^i this 
context, a key issue concerns the build-up of stresses in the crust of these strongly magnetized stars, the eventual 
fracture and the induced dynamics of the coupled magneto-elastic system. These are exciting problem of immediate 
astrophysical interest, and we hope to make progress on them in the near future. 
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Appendix A: Appendix: The non-dimensionalised problem 
1. Dimensionless Background Problem 

The equation of state ([77)) allows us to write the equilibrium equations, ([5]) and in a dimensionless form, c.f., 

P = Pc^^ , (Al) 

P = Pc^^+^ , (A2) 

m = rriQr] , (A-3) 

r = roC, (A4) 

where pc is the density at the centre of the star, and Pc = Kpc is the central pressure. The units of mass become 

rriQ = 47rrgPc: and radius, Tq = ^ ^^ ' where we define b = Pc/Pc- As discussed in Q, 6 is a parameter that may be 
used to gauge the relativistic behaviour of a model. In the limit b —> the model becomes Newtonian. 
Using the dimensionless variables the background equations (O and ([7]) are re- written as [3[ , 

i = et)", (A5) 

where / = 1 — 2(A^ + l)bri/(. Using these variables and considering regularity at the centre, we have the conditions 

For numerical accuracy we change variables one last time. We use x — v/C^ and x = In^ [3|. By using these 
coordinates, we allow for greater accuracy near the centre of the star. The final form of the background equations 
becomes 

X = z?^ - 2x , (A7) 

ij = - y (x + b^^+^) (1 + 6^) , (A8) 

where / = 1 — 2(n + 1)6C^X- Due to difficulties integrating from the actual origin of the star, we start the integration at 
X = —30 which corresponds to C ~ 10^^^ . The integration terminates at C = C/ where ^ drops below a predetermined 
numerical tolerance. With these definitions we define the mass, M, and radius, R, of the star to be Q 



R ^ ^^Ltli^A'/2;,(l-^)/2^^ ^ (AlO) 

where x/ = x(C/k 

We note that Q chose to use b in place of the central density p^ to label a stellar model. This is possible since 
Pc = b^ jK^ where they use {K,N) to parameterize their equation of state. Motivated by the need to determine 

a physical location of the crust, we choose to parameterize our equation of state using {pc,N) using K = bj pj . 
Thus we can specify the central density, as well as the crust-core and crust-ocean transition densities. This choice of 
parameterization does not impact the form of the stellar compactness or the fluid tidal Love number as defined in [3[ . 
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2. Dimensionless Crustal Perturbation Equations 

We perforin the same procedure used in Sec. lA II to reduce our problem to a dimensionless set of equations. First 
we note that our perturbation ^° has units of distance, so we define V -^ rf^V and W -^ r^W. We also note that the 
shear modulus fi has units of pressure, so we define fi — > p^jl. The factor h which should be present is absorbed into 
our definition of jl. 

Using this information in our definition of the traction reveals 



Tr = f Pc.^ 



C^{K - H2) - l{l + l)V - 2CW + (4 - C\)W 



n 



-Afip.rliCV - 2V + e^W) , 



(All) 
(A12) 



where the dots denotes differentiation with respect to C. For notational simplicity we also define SP = TqSP. 
This gives us the dimensionless form of the crustal equations; 



2 X\ ~ C(K - Ho) , , ~ 1(1 + 1) ~ 3Ti 47r 

^ - (c ~ 2 j ^ + ^^^ - «^(^ + ^)^^^ - V^^ - 8A^ (ArTT)5 

C C 4K {N + 1)6 ' 

k = uH + H+ ^(Ci> + 2)V{N + 1)5 - Stt^ , 



T, 



H 



\-v 1 

2 C_ 


T2~'^ 


2 u-\ 

C+ 2 J 


H + 



T2 ~ 2e^C5P + ^{i> + \)H + ^e^[2 - l{l + 1)]V 



1 

8^" 



47r 



C 






-l{l + l)e^ 2(1 -e^) A + 3;> 



e 



e C 



e^{S + l/c4)6P + 16fl 



AN + l)b 



47r 



l-e^ + C(!> + 5A) 



H 



V ,.{N + l)b/:~ 



1>T2 

T 



Here we used p^Tq = {N + l)b/4TT. We also re-write the algebraic relations for Ti (|43p and 6P (|33]) as 



Ti = 



(N + I)b8fl 



4^C' 



4. 3 v(^^+im.+i)R'''+^^'^-^'^^+'^^n'"i'^' ' 



(A13) 

(A14) 
(A15) 

(A16) 



(A17) 



(A18) 
(A19) 



Wne^C^SP ^ C^i>H +[2- 1(1 + l)]e^K + [l{l + l)e^ -2 + C^v'^]H 

+ 8(iV + \)bp[C,'^v^ + [l{l + 1) - 2](1 - e^)]V - Stt{Cv + 2)T2 + 87r(e^ - 3)ri 

Using the same computational variables as in Q for the crust, we eliminate division by the radial coordinate. This 
is more useful for calculations that include the centre of the star; however, it is a simple coordinate transform that 
allows us to join the fluid core to the crust without changing coordinates. 

er^-{K-Ho) 



d^W 



2-^-^\W 



8~^{N + l)be-^^V-^^^±^V ^^' 



An 



8/i {N+l)b 



To 47r 



d.^V = 2V -e^W - ,. ,^^ 

Ajl{N + l)b ' 

d^K = {d:,v)H + d^H + %jl{d^v + 2)V{N + 1)6 - STrTa , 

{dx\ - dxv) 



d.T2^ 
dlH + 



- 1 



{dxv - dxX) 



T2 - 2e^e~^^6P + ^{Oxi^ + 9,A)iJ + 4/ie^[2 - 1(1 + 1)]V—^^ 
Sn [N + 1)6 

dxH + [-1(1 + l)e^ + 2(1 - e^) + (^^.A + 39.i^) - [d^iyf] H 



e^Ti 



(A20) 

(A21) 
(A22) 
(A23) 



e^(3 + l/c2)e-2-JP 



fiV+2)6 



1 - e^ + {dxv + ia^A) - 



d^iy^ 



V + Adxv 



{N + l)b 
47r 



^^{pv) +dxvT2 



(A24) 
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The algebraic relations are also slightly modified to 



IGne^e-^^'SP = <9^i^<9^i7 + [2 - l{l + l)]e^K + [l{l + l)e^ - 2 + id^vf]H 

+ 8{N + l)bfi[{d^iy)^ + [1(1 + 1) - 2](1 - e^)]V - 8Tr{d^iy + 2)7^2 + 87r(e^ - 3)ri . (A26) 

We note that, with our particular formulation it is not obvious that we may consider the Newtonian limit since 
(jA20p and (|A21[) are divergent in that limit. To obtain a non-divergent form we would need to use the equations that 
follow from the divergence of the stress-energy. Since the Newtonian limit is not the focus of our study we do not 
pursue this form of the equations. We also note that the given dimensionless analysis is only valid for a poly tropic 
equation of state. 
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